#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define abs(x) (x >= 0 ? x : -x)
#define max(x, y) (x > y ? x : y)

typedef struct {
    int width;
    int height;
    int stride;

    float * elements;
} Matrix;

int main(int argc, char * argv[])
{
    FILE *ifp; /* Variavel que contem o ponteiro para arquivo de entrada */

    int order = 0; /* Variavel que indica a ordem da matriz */
    int row_test = 0; /* Variavel que indica a linha a ser avaliada */
    int i; /* Variavel contadora */
    int j; /* Variavel contadora */
    int ite_max = 0; /* Variavel que indica o numero maximo de iteracoes */
    int iter = 0; /* Variavel que indica o numero de iteracoes */

    float error = 0.0; /* Variavel que indica erro para a condicao de parada */
    float calc_error = 0.0; /* Variavel do erro calculado */
    float maxl; /* Variavel que indica o maior valor na linha */
    float maxt = 0.0; /* Variavel que indica o maior valor total */
    float res = 0; /* Variavel que contem o resultado final */

    Matrix x0; /* Valores dos x0 */
    Matrix x1; /* Valores dos x1 */
    Matrix mat_a; /* Variavel que contem a matriz A */
    Matrix mat_a_e; /* Variavel que contem a matriz A estrela */
    Matrix mat_b; /* Variavel que contem a matriz B */
    Matrix mat_b_e; /* Variavel que contem a matriz B estrela */

    /* Verifica o uso correto do programa */
    if(argc != 2)
    {
        printf("Usage: ./sequencial [inputfile]\n");
        return 1;
    }

    /* Verifica se foi possivel abrir o arquivo */
    if((ifp = fopen(argv[1], "r")) == NULL)
    {
        printf("Could not open the file %s\n", argv[1]);
    }

    /* Recebe valores de definicao da matriz */
    fscanf(ifp, "%d\n%d\n%f\n%d", &order, &row_test, &error, &ite_max);

    /* Aloca memoria para a matriz A */
    mat_a.elements = (float *) malloc(order * order * sizeof(float));

    /* Aloca memoria para a matriz A estrela */
    mat_a_e.elements = (float *) malloc(order * order * sizeof(float));

    /* Aloca memoria para a matriz B */
    mat_b.elements = (float *) malloc(order * sizeof(float));

    /* Aloca memoria para a matriz B estrela */
    mat_b_e.elements = (float *) malloc(order * sizeof(float));

    /* Aloca memoria para os valores de x0 */
    x0.elements = (float *) malloc(order * sizeof(float));

    /* Aloca memoria para os valores de x1 */
    x1.elements = (float *) malloc(order * sizeof(float));

    /* Recebe os valores do arquivo para a matriz A */
    for(i = 0; i < order; i++)
        for(j = 0; j < order; j++)
            fscanf(ifp, "%f", &mat_a.elements[(i * order) + j]);

    /* Recebe os valores do arquivo para a matriz B */
    for(i = 0; i < order; i++)
        fscanf(ifp, "%f", &mat_b.elements[i]);


    /****************************************************/
    /*************** METODO JACOBI-RICHARDSON ***********/
    /****************************************************/

    /* Calcula o criterio de convergencia e normaliza a matriz A e B */
    for(i = 0; i < order; i++)
    {
        maxl = 0;
        for(j = 0; j < order; j++)
        {
            if(i != j)
            {
                mat_a_e.elements[(i * order) + j] =
                    mat_a.elements[(i * order) + j] /
                    mat_a.elements[(i * order) + i];
                maxl += mat_a_e.elements[(i * order) + j];
            }
            else mat_a_e.elements[(i * order) + j] = 0;
            maxt = max(maxt, maxl);
        }

        mat_b_e.elements[i] = mat_b.elements[i] /
            mat_a.elements[(i * order) + i];
        x0.elements[i] = mat_b_e.elements[i];
    }

    if(maxt > 1)
        printf("Convergence criterion not satisfied.\n");
    else
        printf("Convergence criterion satisfied.\n");

    /* Inicio do algoritmo */
    do {
        iter++; /* Soma uma iteração */

        /*
         * Removido a copia de x1 para x0 para acelerar o processo
         * Agora, se a iteração é ímpar, x1 e x0 representam os seus papeis
         * Quando, a iteração é par, x1 e x0 tem papeis trocados
         */
        if(iter % 2 == 1) /* Iterações ímpares */
        {
            /* Calcula a proxima iteração */
            for(i = 0; i < order; i++)
            {
                x1.elements[i] = mat_b_e.elements[i];
                for(j = 0; j < order; j++)
                    if(i != j)
                        x1.elements[i] -= mat_a_e.elements[(i * order) + j]
                            * x0.elements[j];
            }

            /* Calcula o erro */
            calc_error = (abs(x1.elements[0] - x0.elements[0]) /
                    abs(x0.elements[0]));
            for(i = 1; i < order; i++)
                calc_error = max(calc_error,
                        (abs(x1.elements[i] - x0.elements[i]) /
                         abs(x0.elements[i])));
        }
        else /* Iterações pares */
        {
            /* Calcula a proxima iteração */
            for(i = 0; i < order; i++)
            {
                x0.elements[i] = mat_b_e.elements[i];
                for(j = 0; j < order; j++)
                    if(i != j)
                        x0.elements[i] -= mat_a_e.elements[(i * order) + j]
                            * x1.elements[j];
            }

            /* Calcula o erro */
            calc_error = (abs(x0.elements[0] - x1.elements[0]) /
                    abs(x1.elements[0]));
            for(i = 1; i < order; i++)
                calc_error = max(calc_error,
                        (abs(x0.elements[i] - x1.elements[i]) /
                         abs(x1.elements[i])));
        }

    } while(calc_error > error && iter < ite_max);

    /* O bloco if-else segue a mesma ideia do bloco anterior */
    /* Efetua o calculo do resultado */
    if(iter % 2 == 1)
        for(i = 0; i < order; i++)
            res += mat_a.elements[(row_test * order) + i] * x1.elements[i];
    else
        for(i = 0; i < order; i++)
            res += mat_a.elements[(row_test * order) + i] * x0.elements[i];

    /* Impressao da saida */
    printf("-----------------------------------------\n");
    printf("Iterations: %d\n", iter);
    printf("RowTest: %d => [%f] =? [%f]\n",
            row_test, res, mat_b.elements[row_test]);
    printf("-----------------------------------------\n");

    /* Limpa a memoria alocada e o arquivo aberto */
    free(mat_a.elements);
    free(mat_a_e.elements);
    free(mat_b.elements);
    free(mat_b_e.elements);
    free(x0.elements);
    free(x1.elements);
    fclose(ifp);
    return 0;
}
